home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / math / alged34.zip / ALGEDSRC.ZIP / ALGCALC.C next >
C/C++ Source or Header  |  1996-06-06  |  17KB  |  556 lines
  1. /*--------------------------------------------------------------------
  2.    Alged:  Algebra Editor    henckel@vnet.ibm.com
  3.  
  4.    Copyright (c) 1994 John Henckel
  5.    Permission to use, copy, modify, distribute and sell this software
  6.    and its documentation for any purpose is hereby granted without fee,
  7.    provided that the above copyright notice appear in all copies.
  8. */
  9. #include "alged.h"
  10. /*-----------------------------------------------------------------
  11.    prime factor - replace a number with prime factors
  12. */
  13. void primefact(node *p) {
  14.   int i;
  15.   double v,f;
  16.  
  17.   for (i=0; i<p->nump; ++i)
  18.     primefact(p->parm[i]);
  19.  
  20.   if (p->kind==NUM && !!(v=fabs(p->value)) && whole(v)) {
  21.     if (p->value<0) {
  22.       p->kind = MUL;
  23.       p->nump = 2;
  24.       p->lf = newnum(v);
  25.       p->rt = newnum(-1);
  26.       p = p->lf;
  27.     }
  28.     for (f=2; f<=maxrat && f<v; ++f)
  29.       if (whole(v/f)) {
  30.         v /= f;
  31.         p->kind = MUL;
  32.         p->nump = 2;
  33.         p->lf = newnum(v);
  34.         p->rt = newnum(f);
  35.         p = p->lf;
  36.         f = 1;         /* start loop at the beginning */
  37.       }
  38.   }
  39. }
  40.  
  41. /*--------------------------------------------------------------------
  42.    rational search - this searches for a denominator to a number.
  43. */
  44. int rational_search(node *p) {
  45.   int i,r=0,s=1;
  46.   static double v,n,x,d,n1,d1,x1;
  47.   double err = pow(10,-sigdig);
  48.  
  49.   twirl();
  50.   if (p->kind==NUM) {
  51.     v = fabs(p->value);
  52.     if (p->value<0) s=-1;        /* sign */
  53.     x = modf(v,&n);
  54.     if (x > err) {              /* v is not already an integer */
  55.       x1=5;
  56.       for (d=2; d<=maxrat; ++d) {
  57.         modf(d*v+0.5,&n);            /* n is the round(d*v) */
  58.         x = fabs(n/d - v);           /* x is the error */
  59.         if (x < x1) {
  60.           n1 = n; d1 = d; x1 = x;
  61.         }
  62.       }
  63.       if (x1 < err) {          /* convert p to a ratio of n/d */
  64.         if (n1==0) {
  65.           p->value = 0;
  66.         }
  67.         else {
  68.           p->kind = DIV;
  69.           p->nump = 2;
  70.           p->lf = newnum(n1*s);
  71.           p->rt = newnum(d1);
  72.         }
  73.         ++r;
  74.       }
  75.     }
  76.     else {                /* throw away the very small part */
  77.       p->value = n*s;
  78.       ++r;
  79.     }
  80.   }
  81.   return r;
  82. }
  83.  
  84.  
  85.  
  86. /*-----------------------------------------------------------------
  87.    reduce a/b to lowest terms (with positive denominator)
  88. */
  89. int reduce(double *a,double *b) {
  90.   int r=0;
  91.   double i;
  92.  
  93.   if (*b<0) { *a=-*a; *b=-*b; }
  94.   for (i=2; i<=maxrat && i<=fabs(*a) && i<=*b; ++i)
  95.     if (whole(*a/i) && whole(*b/i)) {
  96.       *a /= i; *b /= i;
  97.       ++r; i=1;
  98.     }
  99.   return r;
  100. }
  101.  
  102. /*---------------------------------------------------------------------------
  103.    ration - this converts all the numbers to ratios of integers if
  104.             possible.
  105. */
  106. int ration(node *p) {
  107.   static double v,u,h,n1,d1,w9,pf,pp,rp;
  108.   static int m,j,k,w,n,f;
  109.   int i,r=0;
  110.   char s[30];
  111.  
  112.   for (i=0; i<p->nump; ++i)
  113.     r+=ration(p->parm[i]);
  114.  
  115.   if (p->kind==NUM && !whole(p->value)) {
  116.     u = fabs(p->value);
  117.     v = frexp(u,&m);             /* u ==> v * 2**m */
  118.     if (m > 0) {
  119.       n = sigdig - m*log10(2);   /* convert m to base 10 logarithm */
  120.       if (n < 1) {
  121.         modf(p->value,&p->value);  /* truncate to integer */
  122.         return r+1;
  123.       }
  124.       v = modf(u,&h);
  125.       m = 0;
  126.     }
  127.     else {                   /* u is a small number */
  128.       h = 0;
  129.       v = modf(u,&h); m=0;  /* suppress tiny fractions */
  130.       n = sigdig;
  131.     }
  132.     sprintf(s,"%1.18f",v);      /* we know 0 < v < 1 */
  133.     /* Look for repeating pattern in s */
  134.     for (f=0; f<n; ++f) {
  135.       k = n-f-1;
  136.       for (w=1; w<k; ++w) {
  137.         for (i=0; i<w; ++i) {
  138.           for (j=f+i+w; j<n && s[j+2]==s[f+i+2]; ) j+=w;
  139.           if (j<n) break;    /* failed */
  140.         }
  141.         if (i==w) break;   /* success */
  142.       }
  143.       if (w<k) break;   /* success */
  144.     }
  145.     if (w<k) {          /* success */
  146.       s[f+w+2] = 0;
  147.       sscanf(s+f+2,"%lf",&rp);     /* repeating part */
  148.       s[f+2] = 0;
  149.       s[1] = '0';
  150.       sscanf(s+1,"%lf",&pp);       /* prefix part */
  151.       w9 = pow(10,w)-1;           /* w nines */
  152.       pf = pow(10,f);
  153.       n1 = (pf*h + pp)*w9 + rp;
  154.       d1 = w9*pf*pow(2,-m);
  155.       if (p->value<0) n1 = -n1;
  156.       reduce(&n1,&d1);
  157.       if (n1==0) p->value = 0;
  158.       else {
  159.         p->kind = DIV;
  160.         p->nump = 2;
  161.         p->lf = newnum(n1);
  162.         p->rt = newnum(d1);
  163.       }
  164.       ++r;
  165.     }
  166.     else             /* try to convert using the search method */
  167.       r += rational_search(p);
  168.   }
  169.   return r;
  170. }
  171.  
  172.  
  173. /*--------------------------------------------------------------------
  174.    Move numbers within clag.
  175. */
  176. int movenums(node *p,int up,int oper) {
  177.   int i,r=0;
  178.   node *a,*b;
  179.  
  180.   for (i=0; i<p->nump; ++i)
  181.     r+=movenums(p->parm[i],up,oper);
  182.   if (p->kind!=oper) return r;
  183.   a = p->lf;
  184.   b = p->rt;
  185.   i = 1;
  186.   if (a->kind!=oper) {
  187.     a=p; i=0;
  188.   }
  189.  
  190.   if (a->parm[i]->kind==NUM) {
  191.     if (p->rt->kind==NUM) {           /* combine nums */
  192.       if (oper==ADD) a->parm[i]->value += p->rt->value;
  193.       else           a->parm[i]->value *= p->rt->value;
  194.       a = p->lf;
  195.       nodecpy(p,a);
  196.       freenode(a);
  197.       freenode(b);
  198.       ++r;
  199.     }
  200.     else if (up) {                         /* switch */
  201.       p->rt = a->parm[i];
  202.       a->parm[i] = b;
  203.       ++r;
  204.     }
  205.   }
  206.   else if (b->kind==NUM && !up) {   /* switch */
  207.     p->rt = a->parm[i];
  208.     a->parm[i] = b;
  209.     ++r;
  210.   }
  211.   return r;
  212. }
  213. /*-----------------------------------------------------------------
  214.    is complex?  returns >0 and a and b if p is a complex (1=real)
  215.    a is the imaginary part.  sorry, this is against convention!
  216. */
  217. int iscomplex(node *p,double *a,double *b) {
  218.  
  219.   if (aop(p->kind) && p->rt->kind==NUM &&              /* ai+b */
  220.       p->lf->kind==MUL && p->lf->rt->kind==VAR &&
  221.       p->lf->lf->kind==NUM && !strcmp(p->lf->rt->name,"i")) {
  222.     *a = p->lf->lf->value;
  223.     *b = p->rt->value;
  224.     if (p->kind==SUB) *b = -*b;
  225.     return 5;
  226.   }
  227.   if (p->kind==MUL && p->rt->kind==VAR &&              /* ai */
  228.       p->lf->kind==NUM && !strcmp(p->rt->name,"i")) {
  229.     *a = p->lf->value;
  230.     *b = 0;
  231.     return 3;
  232.   }
  233.   if (aop(p->kind) && p->rt->kind==NUM &&              /* ia+b */
  234.       p->lf->kind==MUL && p->lf->lf->kind==VAR &&
  235.       p->lf->rt->kind==NUM && !strcmp(p->lf->lf->name,"i")) {
  236.     *a = p->lf->rt->value;
  237.     *b = p->rt->value;
  238.     if (p->kind==SUB) *b = -*b;
  239.     return 5;
  240.   }
  241.   if (p->kind==MUL && p->lf->kind==VAR &&              /* ia */
  242.       p->rt->kind==NUM && !strcmp(p->lf->name,"i")) {
  243.     *a = p->rt->value;
  244.     *b = 0;
  245.     return 3;
  246.   }
  247.   if (aop(p->kind) && p->lf->kind==VAR &&            /* i+b */
  248.       p->rt->kind==NUM && !strcmp(p->lf->name,"i")) {
  249.     *a = 1;
  250.     *b = p->rt->value;
  251.     if (p->kind==SUB) *b = -*b;
  252.     return 4;
  253.   }
  254.   if (p->kind==VAR && !strcmp(p->name,"i")) {          /* i */
  255.     *a = 1;
  256.     *b = 0;
  257.     return 2;
  258.   }
  259.   if (p->kind==NUM) {                              /* b */
  260.     *a = 0;
  261.     *b = p->value;
  262.     return 1;
  263.   }
  264.   return 0;
  265. }
  266.  
  267. /*-----------------------------------------------------------------
  268.    complex base raised to a power - convert to r^x * cos tx + i sin tx
  269.    b is the imaginary part.  p->rt is REUSED and p->lf is freed.
  270.     -----------------------------------------------------------------
  271.          Note: I am not using this function because (1) it is
  272.          questionable whether the cos+isin expansion is "simpler" than
  273.          exp(i*x), and (2) because it was happening needlessly when the
  274.          complex base was not completely sorted and simplified.
  275.  
  276. void complexpow(node*p,double b,double a) {
  277.   node *c,*s;
  278.   double t;
  279.  
  280.   t = atan2(b,a);         //  -pi < t <= pi
  281.   c = newnode();
  282.   c->kind=FUN;
  283.   strcpy(c->name,"cos");
  284.   c->nump=1;
  285.   c->lf = cons(deepcopy(p->rt),MUL,newnum(t));
  286.  
  287.   s = newnode();
  288.   s->kind=FUN;
  289.   strcpy(s->name,"sin");
  290.   s->nump=1;
  291.   s->lf = deepcopy(c->lf);
  292.  
  293.   freetree(p->lf);
  294.   p->kind = MUL;
  295.   p->lf = cons(newnum(hypot(a,b)),EXP,p->rt);
  296.   p->rt = cons(c,ADD,cons(newvar("i"),MUL,s));
  297. }
  298. ****************/
  299.  
  300. /*--------------------------------------------------------------------
  301.    calcnode
  302.    calculate as many numeric results as possible
  303. */
  304. int calcnode(node *p,int keeprat) {
  305.   int i,r=0;
  306.   node *a,*b;
  307.   static double a1,b1,a2,b2,r1,t1,r2,t2;
  308.  
  309.   for (i=0; i<p->nump; ++i)
  310.     r+=calcnode(p->parm[i],keeprat);
  311.  
  312.   a = p->lf;
  313.   b = p->rt;
  314.   switch (p->kind) {
  315.   case ADD:
  316.     if (a->kind==NUM && a->value==0) {          /* 0+x = x */
  317.       nodecpy(p,b);
  318.       freenode(a); freenode(b);
  319.     }
  320.     else if (b->kind==NUM && b->value==0) {     /* x+0 = x */
  321.       nodecpy(p,a);
  322.       freenode(a); freenode(b);
  323.     }
  324.     else if (a->kind==NUM && b->kind==NUM) {    /* n+n = n */
  325.       p->kind = NUM;
  326.       p->nump = 0;
  327.       p->value = a->value + b->value;
  328.       freenode(a); freenode(b);
  329.     }
  330.     else if (b->kind==MUL &&                   /* a+-1*b = a-b */
  331.              b->lf->kind==NUM &&
  332.              b->lf->value<0) {
  333.       p->kind = SUB;
  334.       b->lf->value = -b->lf->value;
  335.     }
  336.     else if (b->kind==MUL &&                   /* a+b*(-1) = a-b */
  337.              b->rt->kind==NUM &&
  338.              b->rt->value<0) {
  339.       p->kind = SUB;
  340.       b->rt->value = -b->rt->value;
  341.     }
  342.     else break; ++r;
  343.     break;
  344.   case SUB:
  345.     if (a->kind==NUM && a->value==0) {          /* 0-x = -x */
  346.       p->kind = MUL;
  347.       a->value = -1;
  348.     }
  349.     else if (a->kind==NUM && b->kind==NUM) {    /* n-n = n */
  350.       p->kind = NUM;
  351.       p->nump = 0;
  352.       p->value = a->value - b->value;
  353.       freenode(a); freenode(b);
  354.     }
  355.     else if (b->kind==NUM && b->value==0) {     /* x-0 = x */
  356.       nodecpy(p,a);
  357.       freenode(a); freenode(b);
  358.     }
  359.     else break; ++r;
  360.     break;
  361.   case MUL:
  362.     if (a->kind==NUM && b->kind==NUM) {         /* n*n = n */
  363.       p->kind = NUM;
  364.       p->nump = 0;
  365.       p->value = a->value * b->value;
  366.       freenode(a); freenode(b);
  367.     }
  368.     else if (a->kind==NUM && a->value==0) {     /* 0*x = 0 */
  369.       p->kind = NUM;
  370.       p->nump = 0;
  371.       p->value = 0;
  372.       freenode(a); freetree(b);
  373.     }
  374.     else if (b->kind==NUM && b->value==0) {     /* x*0 = 0 */
  375.       p->kind = NUM;
  376.       p->nump = 0;
  377.       p->value = 0;
  378.       freetree(a); freenode(b);
  379.     }
  380.     else if (a->kind==NUM && a->value==1) {     /* 1*x = x */
  381.       nodecpy(p,b);
  382.       freenode(a); freenode(b);
  383.     }
  384.     else if (b->kind==NUM && b->value==1) {     /* x*1 = x */
  385.       nodecpy(p,a);
  386.       freenode(a); freenode(b);
  387.     }
  388.     else if (iscomplex(a,&a1,&b1) && iscomplex(b,&a2,&b2) && // (ai+b)*(ci+d)
  389.             !!(t1 = b1*b2 - a1*a2)) {
  390.       freetree(a);
  391.       freetree(b);
  392.       p->rt = newnum(t1);
  393.       p->kind = ADD;
  394.       p->lf = cons(newnum(b1*a2 + b2*a1),MUL,newvar("i"));
  395.     }
  396.     else break; ++r;
  397.     break;
  398.   case DIV:
  399.     if (b->kind==NUM && b->value==0) {          /* x/0 = BAD */
  400.       p->kind = BAD;
  401.       p->nump = 0;
  402.       p->value = 0;
  403.       freetree(a); freenode(b);
  404.     }
  405.     else if (a->kind==NUM && a->value==0) {     /* 0/x = 0 */
  406.       p->kind = NUM;
  407.       p->nump = 0;
  408.       p->value = 0;
  409.       freenode(a); freetree(b);
  410.     }
  411.     else if (b->kind==NUM && b->value==1) {     /* x/1 = x */
  412.       nodecpy(p,a);
  413.       freenode(a); freenode(b);
  414.     }
  415.     else if (a->kind==NUM && b->kind==NUM) {    /* n/n = n */
  416.       if (keeprat &&
  417.           whole(a->value) &&       /* if both int, then reduce */
  418.           whole(b->value)) {
  419.         if (!reduce(&a->value,&b->value)) break;
  420.       }
  421.       else {
  422.         p->kind = NUM;
  423.         p->nump = 0;
  424.         p->value = a->value / b->value;
  425.         freenode(a); freenode(b);
  426.       }
  427.     }
  428.     else if (iscomplex(a,&a1,&b1) && iscomplex(b,&a2,&b2)) { // (ai+b)/(ci+d)
  429.       r2 = a2*a2 + b2*b2;
  430.       if (!r2) break;              // divide by zero, no change
  431.       movenode(p,cons(cons(newnum((a1*b2 - a2*b1)/r2),MUL,newvar("i")),
  432.                       ADD,newnum((a1*a2 + b1*b2)/r2)));
  433.     }
  434.     /*-----------------------------------------------------------------
  435.        the following is called "stretch rule" to accomidate numbers
  436.        separated by the bisection function.
  437.     */
  438.     else if (
  439.         (a->kind==NUM || a->kind==MUL && a->rt->kind==NUM && !!(a=a->rt)) &&
  440.         (b->kind==NUM || b->kind==MUL && b->rt->kind==NUM && !!(b=b->rt))) {
  441.       if (keeprat &&
  442.           whole(a->value) &&       /* if both int, then reduce */
  443.           whole(b->value)) {
  444.         if (!reduce(&a->value,&b->value)) break;
  445.       }
  446.       else {
  447.         a->value = a->value / b->value;
  448.         b->value = 1.0;
  449.       }
  450.     }
  451.     else if (b->kind==NUM ||
  452.              b->kind==MUL && b->rt->kind==NUM && !!(b=b->rt)) {
  453.       if (keeprat && whole(b->value)) {
  454.         if (b->value >= 0) break;
  455.         p->lf = cons(p->lf,MUL,newnum(-1));
  456.         b->value = -b->value;
  457.       }
  458.       else {
  459.         p->lf = cons(p->lf,MUL,newnum(1.0/b->value));
  460.         b->value = 1.0;
  461.       }
  462.     }
  463.     else break; ++r;
  464.     break;
  465.   case EXP:
  466.     if (b->kind==NUM && b->value==0) {          /* x^0 = 1 */
  467.       p->kind = NUM;
  468.       p->nump = 0;
  469.       p->value = 1;
  470.       freetree(a); freenode(b);
  471.     }
  472.     else if (a->kind==NUM && a->value==1) {     /* 1^x = 1 */
  473.       p->kind = NUM;
  474.       p->nump = 0;
  475.       p->value = 1;
  476.       freenode(a); freetree(b);
  477.     }
  478.     else if (b->kind==NUM && b->value==1) {     /* x^1 = x */
  479.       nodecpy(p,a);
  480.       freenode(a); freenode(b);
  481.     }
  482.     else if (iscomplex(a,&a1,&b1) &&
  483.              iscomplex(b,&a2,&b2)) {           // (ai+b)^(ci+d)
  484.       if (keeprat && whole(a1) && whole(b1) &&
  485.                      (a2 || !whole(b2) || b2<0)) break;  // no change
  486.  
  487.       //  First we check for some simple cases because the
  488.       //  general method has problems with round-off errors.
  489.       if (!a2 && !a1 && b1>0)           // positive real base
  490.         movenode(p,newnum(pow(b1,b2)));
  491.       else if (!a2 && !a1 && whole(2*b2))  // (-x)^(n/2) = i^n * x^(n/2)
  492.         movenode(p,cons(cons(newvar("i"),EXP,newnum(2*b2)),
  493.                         MUL,newnum(pow(-b1,b2))));
  494.       else if (!a2 && !b1 && whole(b2)) {  // (ai)^n = a^n * {1,i,-1,-i}
  495.         i = b2;
  496.         if (i&1)
  497.           movenode(p,cons(newvar("i"),MUL,newnum(pow(a1,b2)*(1-(i&2)))));
  498.         else
  499.           movenode(p,newnum(pow(a1,b2)*(1-(i&2))));
  500.       }
  501.       else {              // general complex exponentiation
  502.         r1 = hypot(a1,b1);
  503.         r2 = hypot(a2,b2);
  504.         if (r2==0)
  505.           movenode(p,newnum(1));
  506.         else if (r1==0)
  507.           movenode(p,newnum(0));     // need this to avoid domain errors
  508.         else {
  509.           t1 = atan2(a1,b1);
  510.           r2 = pow(r1,b2)*exp(-t1*a2);
  511.           t2 = t1*b2 + a2*log(r1);
  512.           movenode(p,cons(cons(newnum(r2*sin(t2)),MUL,newvar("i")),
  513.                       ADD,newnum(r2*cos(t2))));
  514.         }
  515.       }
  516.     }
  517.     else break; ++r;
  518.     break;
  519.   case FUN:
  520.     if (p->nump==1 && a->kind==NUM) {
  521.       if (!strcmp(p->name,"sin")) p->value=sin(a->value);
  522.       else if (!strcmp(p->name,"cos")) p->value=cos(a->value);
  523.       else if (!strcmp(p->name,"tan")) p->value=tan(a->value);
  524.       else if (!strcmp(p->name,"acos")) p->value=acos(a->value);
  525.       else if (!strcmp(p->name,"asin")) p->value=asin(a->value);
  526.       else if (!strcmp(p->name,"atan")) p->value=atan(a->value);
  527.       else if (!strcmp(p->name,"cosh")) p->value=cosh(a->value);
  528.       else if (!strcmp(p->name,"sinh")) p->value=sinh(a->value);
  529.       else if (!strcmp(p->name,"tanh")) p->value=tanh(a->value);
  530.       else if (!strcmp(p->name,"ln")) p->value=log(a->value);
  531.       else if (!strcmp(p->name,"log")) p->value=log10(a->value);
  532.       else if (!strcmp(p->name,"abs")) p->value=fabs(a->value);
  533.       else if (!strcmp(p->name,"rand")) p->value=a->value*rand()/RAND_MAX;
  534.       else if (!strcmp(p->name,"sign")) p->value=a->value<0?-1:a->value>0?1:0;
  535.       else break;
  536.       p->kind = NUM;
  537.       p->nump = 0;
  538.       freenode(a);
  539.     }
  540.     if (p->nump==2 && a->kind==NUM && b->kind==NUM) {
  541.       if (!strcmp(p->name,"max")) p->value=max(a->value,b->value);
  542.       else if (!strcmp(p->name,"min")) p->value=min(a->value,b->value);
  543.       else if (!strcmp(p->name,"atan2")) p->value=atan2(a->value,b->value);
  544.       else if (!strcmp(p->name,"mod")) p->value=fmod(a->value,b->value);
  545.       else if (!strcmp(p->name,"r")) p->value=hypot(a->value,b->value);
  546.       else break;
  547.       p->kind = NUM;
  548.       p->nump = 0;
  549.       freenode(a); freenode(b);
  550.     }
  551.     break;
  552.   }
  553.   return r;
  554. }
  555.  
  556.